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CHAPTER 2 


RADIATIVE TRANSFER THROUGH THE ATMOSPHERE 

In this chapter the numeric solution to the transfer of visible 
and near infrared energy within the atmosphere is discussed. Models of 
the atmosphere are included. Particular emphasis is given to Rayleigh and 
Mie theories of scattering, as well as absorption due to atmospheric 
ozone and water. 

The fundamental mathematics for the theory of radiative transfer 
were developed by Chandrasekhar (1950). He was first to formalize the 
problem of radiative transfer in a solar illuminated plane-parallel 
atmosphere, and present a solution in the form of a set of nonlinear 
equations. He accounted for polarization by adopting the four Stokes 
parameters to characterize the field. However, even with this framework 
the solution to these equations remained unsolved for several years. This 
can be attributed to the fact that the radiance, for a given altitude 
within the atmosphere and directed in a given direction, is expressed in 
terms of the radiances incoming from all directions, for that altitude. 

In practice, the closed form solution, which expresses the field in terms 
of the known boundary values, cannot be written. The number of equations 
involved would be overwhelming. Approximations, such as that of single 
scattering, are often made. This is acceptable provided high accuracy is 
not a criteria. The LOWTRAN 6 code (Kneizys and associates, 1983), 
available through the Air Force Geophysics Laboratory, is one such 
approximate code that is based on the assumption of single scattering. 



Using the Gsuss-Seidel Iterative technique, Herman (1963; Herman 
and Browning, 1965; Herman, Browning, and Curran, 1971) developed the 
software required to solve the Chandrasekhar equations. In his method an 
Initial guess Is made of the field present after passing through one 
layer, A solution for successive layers Is made using quantities that 
have been calculated In the previous layer. At the ground a reflectance 
model, usually lambertlan, Is used to compute the upwelling radiance* 
Radiances are then traced moving back up to the top of the atmosphere. 
Once the radiances at all atmospheric levels have been solved, the 
process Is repeated utilizing updated values of the assumed radiances. 
All unknowns are changed from their previously calculated values, as the 
value of the Initial unknowns are changed. After several iterations the 
unknowns converge to a unique solution. This solution is exact In the 
sense that all orders of multiple scattering are accounted for. 
Computational accuracy is thought to be limited only by the atmospheric 
models and input parameters required to run the code. 

We use this radiative transfer program for use in the absolute 
radiometric calibration of Landsat’s Thematic Mapper. In the following it 
is referred to as the Herman Code. 

The Equation of Transfer 

The attenuation of radiation through some distance ds can be 
described by the equation 

dL^ = “k-jx p Lx ds * (2.1) 

Here kxx is the total mass extinction coefficient (in units of area per 
mass, such as cra z /gra), p is the density of the medium (mass/volume, or 

gm/cra 4 ), and L is the radiance (W/cm J sr pm) at point s within the medium. 
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The distance ds is a positive quantity, irrespective of coordinate system.. 
Due to this attenuation, and due to scattering into the beam, radiance 
varies with distance a. Extinction and density may also vary spatially. 
In addition, the radiance and mass extinction coefficient vary with 
wavelength. It is common to identify spectral concentrations with a 
subscript X (such as L^) and spectral functions as (x), but for simplicity 
these notations will not be used further. 

The mass extinction coefficient is composed of a scattering term, 
k s » and an absorption term, k a . Thus, 

kx “ k e + k a . (2.2) 

A related parameter is the volume extinction coefficient, kf, which 

has units of inverse length. Usually one prefers to describe the 
variation of extinction within the atmosphere in terms of the particle or 
molecular density. Thus, the radiative transfer equations most often use 
the mass extinction coefficient, as opposed to the volume extinction 
coefficient. The former is usually assumed constant throughout the 
atmosphere. It will only vary spatially with altitude if effects such as 
pressure broadening, variations of aerosol refractive index, or variations 
in aerosol radial size distribution occur. Conversely, the parameter 
varies dramatically with altitude due to its proportionality to density. 

To describe the distribution of radiance, normalized to the 
incoming irradiance, that is scattered from a beam, the phase function 
P(e) is introduced. Here e is the angle between the incident and 
scattered beams. Equivalently, P(e» <$>5 0* describes that radiance 
scattered from a differential solid angle centered about (0,<{>) into a 
differential solid angle about C 0 * , <t> * ) * In this chapter, the first angles 
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within parenthesis will be those of the incident beam; the angles which 
follow the semi-colon are those of the scattered beam. When P(0) is used 
within an integrand, the integration will be with respect to the primed 
angles. For example, by integrating the phase function over all outgoing 
angles, dio 1 , the total energy lost from a beam through scattering can be 
computed. This is given by 

AL g (0,4>) =* -k 0 p L( Q , <t>) ds (2.3) 

2 Trr 7T 

« -kf p ds P( 6 , 4>; 0' ,$') L(0, $) sin 0' d 0 ' d <t>' . 

Jo lo 

L(0,iji) is ,:ot a function of the scattered angles, and may be placed 
outside the integral. Using the above equality, the identity 

/ 2 TTf IT k 

P(e,4); 0* ,4>*) sin e' d0' d<j>' = ^ = u ( (2.4) 

Jo Jo T 

is made. The ratio kg/kj is known as the single scattering albedo, w 0 . 

It is that fraction of the total attenuation due to scattering for a 
single collision, and is equal to the integral of the phase function over 
the scattered angles. A conservative scattering atmosphere is one in 
which w„=l. 

The energy balance equation, which summarizes the sources and 
sinks acting on a beam, can be written 

r 

dL( Q , <}>) = k^ p ds P(0',* , }9,+) L( 0 * , 4> ’ ) dw' (2.5) 

■ 4ir 

+ kT P ds P(0o,<t>o;6, 40 E + ep ds - kT p L(9,4>) ds. 

The first terra is that energy scattered into dm due to incoming fields 

from all directions dm', where du'^sine 1 d 0 ' d<j>'. Note how this integral 
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over Che phase function differs from before. The energy into L(0,<}>) from 
L< 0 ’ , ’ ) is computed here, as opposed to the energy out of L(0,<{>) into all 
LC 0 1 » ' ) as in Equation (2.3), The second term accounts for single 
scattering out of the solar beam. The incident solar beam has an 
irradiance E at distance s and propagates along (0 # ,4> o ) where o 0 >90°, 
0o a 0z + 9O°, and 0 Z is the solar zenith angle. The third terra accounts for 
emission within the atmosphere. The spectral parameter e^ (here denoted 
only as e) is the emitted spectral radiant flux propagating in an 
infiniti8imal cone containing that direction of propagation, divided by the 
solid angle of the cone, and normalized with respect to the density of 
the medium. The final term represents Chat energy lost due to scattering 
and absorption processes. 

The sources can be readily grouped together by introducing the 
source function 


J(0,$) 


P(0',$'; 8, <j>) L(0',<j>') dm' + Pq( Bo » 4>o> 9> E + e/^T* (2.6) 

4tt 


This is the radiance added to the incident beam from a unit mass of the 
medium and for a mass extinction coefficient of unity. Throughout the 
visible and near infrared regions of the spectrum, emission is considered 
to be negligible. For our application, therefore, the source function 
will only have contributions from scattering. After dividing both sides 
by (-kx p ds), and introducing the source function, Equation (2.5) becomes 


dL(8.6) 
-kx p ds 


— L(o,$) — J( Q , <{>) • 


(2.7) 



The derivative with respect to s is now expanded in terms of 
derivatives with respect to the x, y, and z axes. Here a cartesian 
coordinate is defined such that the z-axis is directed upward and the x- 
axis Is directed such that the sun falls within the x-z plane. In 
addition, the zenith angle 0 is defined with respect to an outward normal 
directed along the z-axis. A beam propagating along 0*0° is propagating 
out towards space; a beam directed into the sun will have an azimuth 
angle of $=0°. This coordinate system is depicted in Figure 2.1. 

Several assumptions are placed on tne atmosphere to be modeled. 
The atmosphere is assumed (a)to be in steady-state (no variations with 
time), and (b)horizontally homogeneous, which implies a flat earth. (As 
the following equations use the steady-state assumption, they cannot be 
used to describe the propagation of a pulsed lidar beam. Here significant 
changes occur within the atmosphere as the beam travels.) These 
assumptions imply that there will be no variations in the field or source 
function along a horizontal plane, and there will be no variations with 
time. Thus if the generalized function f here represents either L or J, 
the derivatives df/dx, df/dy, and df/dt will be zero and df/ds « 
(df/dz)(dz/ds) = (df/dz) cos©. 

A few definitions may be conveniently introduced here. First let 
p=|cos0|. With the sun at a solar zenith angle of $ Zt rays propagating 
downward from the sun are associated with -p 0 =! cos( 0 Z +9O°)| • While p 
itself is always positive, the angle -p will be associated with downward 
propagating beams and p will be associated with upward directed beams. 
Secondly, let the optical depth at any height z within the atmosphere be 
defined as 
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Figure 2.1 Coordinate system used to represent beam directionality. 
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Figure 2.2 Layer nomenclature for beams propagating 
through the atmosphere. 
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k p dz. 
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The optical depth between two altitudes la given as 


At(zj ,Zj) 


’Zl 

k p dz 
Zl 


(2.9) 


- t(z ,) - t(z 2 ). 

As t is a positive quantity which monotonically decreases with increas 
altitude z, the integration within Equation (2.9) will always be set up 
such that z 2 >Zj and hence At(zj,z 2 ) will be positive. If z is at ground 
elevation, and k=k? the optical depth defined by (2.8) is T exC , the 
extinction optical thickness, or total optical depth of the atmosphere. 

The attenuation of radiance can now be written in terms of 
optical depth 

dL(8,<j)) 3 -kx p L(0,$) dz/cos6 3 dt L(0,<j>)/cos0 . (2.10) 

Here the substitution dT“-kT P ds and ds-dz/ccsO have been made within 
Equation (2.1). Note ds=dz/cos0 is again a positive quantity. For 
downward directed beams both dz and coso are negative; for upward 
directed beams both are positive. This is necessary to assure that dL, as 
given by Equation (2.10), is always negative. (The beam is attenuated by 
this term.) Two separate equations are, however, required to express ds 
as a function of p, as p in itself carries no sign. For downward directed 
beams ds=~dz/p, while for upward directed beams ds-dz/p. 

The energy balance equation, Equation (2.7), may likewise be 
rewritten: 
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downward 


upward 


-p dL(T,-p»4>)/dT * L( t ,- p, $) - J(t,-p,$) (2.11a) 


U dL(T,p,<t»)/dT « L(t,u,4.; - J(t,p,$) , (2.11b) 

For the time being it will be convenient to write separate equations for 
downward and upward directed beams. Note the radiance and source terms 
are a function of both altitude z (hence a function of t) and direction 
(p,<}>)« The parameters within parenthesis serve as a reminder of this 
dependence. The equation is a nonlinear, first order differential 
equation, subject tfj the following boundary value conditions. The diffuse 
radiance incident at the top of the atmosphere is zero, there are no 
contributions to radiance from below the earth's surface, and the exo- 
atmospberic solar lrradlance is known. That is, 

L(0,- u ,$) - 0 (2.12) 

Tex t » M * a 0 

" known 


The parameter r e xt *8 the total optical depth at the earth's surface. 

Multiplying both sides of Equation (2.11a) by exp(t/p) and both 
sides of Equation (2.11b) by exp(-x/p) we obtain 
for downward propagation 

-p e T /y dL(f,-p, <J>)/dx “ L(x,-*p,<t>) e 1 /* 1 3 -p — — (2.13a) 

■ ~J(t> _ p» < t>) ef/p 

and for upward propagation 

p e-T/p dL(x,p,4>)/dT - L(x,p,4>) e-t/p = p (2.13b) 

a x 

3 -J(x,p,4>) 
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Consider a ray as It traverses the layer structure shown in Figure 2.2. 
The top of the layer is denoted by x n on^ the bottom by x n +i. These 
layers are also denoted with xj (initial) and Tf (final), where x^ can be 
either x n or T n +,, depending on the direction of propagation. 
Intermediate altitudes are identified by some x'. The radiance after 
propagation is determined by integrating Equation (2.13) between the 
initial and final t values. Thus, 
downward 

-ji [L(T tJ+l ,- l j,<j>) etn+i/l 1 - L(x n ,-p, $) e T n/^] (2.14a) 


( xn+i 


J(x',-M) e T ’/>i dx' 


upward 


V [L(t„ ,!!,♦) e~ T n/» - L(x n+l> p,*) e- T n+i/b] 


(2.14b) 


“ - J(x',p,iji) dx' 

T n+i 

Dividing Equation (2.14a) by (-P exp(x n+l /p)) and Equation (2.14b) by (M 

exp(-x n /p)), the above can be rewritten: 

downward 

L(T n+t ,-H, ♦) ■ L(x n ,-u,<|>) e"( T n+i“ T n)/w (2.15a) 


+ 


T n+i 

j( x f ,y , +) e -, ( T n+l“ T, )/l 1 dx'/n 
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upward 


L(T n ,M) - L(T n+1 ,u»*) e-(Tn+i"Tn)/p 


(2.15b) 


T n+i 


j(x ' ,| 1 ,$) e -(T»“T n )/p dt'/jl . 


; T n 


These equations con be combined if and xf are introduced. 


L(xf,±y,<|>) - L(tj,±M) e“ AT /l J + 


T n+i 


J(x',±y,$) dx’/y (2.16a) 


where 


d(x',±y, 4>) 


F(ji , ,4>';±y, +) L(y* , # ) d(-y') d<j> + P(y 0 ,4> 0 ; 0, ♦) E . (2.16b) 


1 4 TT 

Here AT«T n +,“T n and At ,3 » |t t f | • The solid angle du'-sinO* d 6 1 d$' ° 
d(cos6')d^' has been written in terras of y, or dto' “d(~y') d$'« Each of the 
above equations state that the radiance after passing through a layer can 
be expressed as the initial radiance attenuated by exp(-Ax/y), plus a 
contribution from the source function. The source function adds a 
contribution at each altitude x*, but is attenuated due to the Ax' between 
x' and the final layer. 


Numeric Solution 

To evaluate the above radiative transfer equation, the integrals 
within Equations (2.16a) and (2.16b) are replaced with an equivalent sum 
of integrals with smaller differences between the limits of integration. 
The new limits are defined such that the parameters J, P, and L can be 
approximated as constants within the Ax, A0, and A$ intervals. They are 
put outside the integrals, an evaluation is made, and a solution is 
obtained. 
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INTRODUCTION 


This is the eighth quarterly report on Contract NAS5-27382 entitled 
"Spec troradio me trie Calibration of the Thematic Mapper and the 
Multispectra 1 Scanner System". 

On 26 - 28 October 1984, we made a successful trip to White Sands, 
New Mexico, and took a variety of radiometric measurements on the 
morning of the Landsat 5 Thematic Mapper overpass. The sky that morning 
was cloud free, our sites were dry although many other areas at White 
Sands were covered by several centimeters of water due to the unusually 
rainy Fall. 

Barnes multiband radiometer data were collected for a 4 X 4 pixel 
area and two fractional pixel areas of slightly higher and lower 
reflectances than the larger area. Helicopter color photography was 
obtained of all the ground areas. This photography will allow us to 
make a detailed reflectance map of the 4X4 pixel area and we will be 
able to register it to the TM imagery to an accuracy of better than half 
a pixel. Spectropolarime ter data were also collected of the 4X4 pixel 
area from the helicopter. In addition, ground based solar radiometer 
data were collected to provide spectral extinction optical thickness 
values. The uncorrected Thematic Mapper image, in CCT form, of White 
Sands for that date is expected in a few days. The completion of the 
calibration of the TM for that date awaits receipt of the CCT. A 
description of the data reduction and calibration will be given in the 
next quarterly report. 

The remainder of this report consists of a description of the 
ladiative transfer theory used in the development of the Herman code 
which we use In predicting the TM entrance pupil spectra? radiances from 
the ground based measurements just mentioned. The theory is then 
fundamental to the measurement program we are conducting at White 
Sands. It has been written up by C.J. Kastner and is also to be included 
as a chapter in her Ph.D. dissertation. 



In evaluating the integral over optical depth, Equation (2.16a), 
the radiative transfer between layers x n and T n + 2 is considered. This 
allows the average value of the source term to be taken as that at the 
midpoint of the interval, namely that at T n + t . After factorin out this 


f Tn+2 

constant value and evaluating J( x n+1 ,±u, $) e -Ax'/u dx'/y, Equation 

* T n 

(2,16a) becomes 

L( Tf ,±y, «t>) “ L( ,±W, e -*Ax/ii + J( T n+i*±U, <f>) (l-e” z At / u ) (2.17) 
where the interval Ax is still defined as that between x n and T n+i> and 
Tf and tj are separated by a 2Ax thickness. 

In turn, J(T n -j. t ,±y t <ji) is evaluated by replacing the integrals 
within (2.16b) with sums over the finite differences Ap and A$. That is, 

(2.18) 


J ^ T n-. u i » ± P»4>) 


2 xr/A 4> 

1 . 


k=l 


ff/A9 

£ PCy , j,« , j c ; ± y, <t») L(x tl+l , (1 'j,<). , k ) (-Ap')j 
j=i 


(A^’)k « 


As an example, let A0=1O° and A 4>=30 0 (actually their radian equivalents). 


Then , 

j - 1 18 (2.19) 

0'j - j A0 - A0/2 » [5°, 15°, .... 1 7 5 0 ] 

jj'j = COS 0'j 

(“Ap')j = cos( (j— 1)A0) - cos(jA0) 

= [cosO°-coslO # , ..., cosl70°-cosl80° ] 

k a 1, 12 

({('k = k A 4> - A$/2 = [15°, 45", ..., 345 s ] 

(a $ ' ) k = A4> n 30° * tt rad/180 0 

Note F and L have been taken out of the integral over the finite limits 
A6 and A$, and replaced with their values at the midpoints of these finite 
differences. 
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At the beginning of each iteration through the atmsophere, the 
radiances at level Ti are required. This is achieved by considering the 
transfer of radiation through only a single At layer, 

2ir/A<j> w / A 9 

Ut,,-^*) “ (1-e-At/ji) [ £ £ (2.20) 

k-1 j-l 

(-Ap')j A<j> + PC~Mo » 4>) E»] . 

On the first pass all upwelling radiances L(0 * w'j , 4» ' k) » or that 
energy being reflected out of the atmosphere and into space, are assumed 
zero. On successive passes, those values computed in previous iterations 
are assumed. At all times the downwelling radiances L(0,h'j ,$'k), at the 
top of the atmosphere, are assumed zero. This is a statement that the 
only energy entering the atmosphere is from the solar irradlance E 0 . 
Similarly, at the bottom of the atmosphere the radiances at Text-i ar * 
computed from those at a single At layer below. At x e xt fc he radiances 
are those reflected from the surface. They are computed by multiplying 
the sum of the diffuse and direct downwelling irradiances by p/ir. 

In choosing a numeric value for the layer thickness, At, Herman 
(1963) used a statistical analysis to compute the probability that 
scattering within a layer would be due to single scattering alone. The 
At interval must be small enough to neglect variations of the source 
term, which is equivalent to requiring that L and E remain approximately 
constant over the interval. This is likely if a photon has a small 
probability of undergoing two or more scattering events. Conversely, At 
must not be so small as to make the computation time excessive. A value 
of AT a 0.02 was chosen. Here, approximately 96% of the scattered 
radiation is associated with a single collision. Since the effective depth 
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of the atmosphere is dtp, a greater percentage of multiple scattering 
occurs at larger zenith angles. As 9 approaches 90®, this error builds up 
rapidly. Calculations down to 85® can, however, be made without 
introducing any serious errors. 

Polarization 

The Herman Code which we have used to date does not account for 
changes in polarization as a ray propagates through the atmosphere. The 
code can, however, be easily modified to do bo. Preliminary studies have 
indicated that the nonpolarization code is accurate enough, given the 
atmospheric conditions we have encountered at White Sands to date, for 
our calibration work* For this reason the studies within this 

dissertation have been made using the original Herman code, which is both 
easier and faster to run. For completeness, the theory behind the 
polarization code is discussed here. 

To be as accurate as possible, the radiative transfer equation 
must describe the state of polarization of a scattered field, as this 
field generally has undergone a change in polarization compared to that 
of the incident field. To describe this state, the amplitude of the 
electric field components along two orthogonal directions, as well as the 
phase difference between these components, are required. For example, 
let E2 and E r be the parallel, and perpendicular, components, defined with 
respect to a reference plane. This reference plane is chosen as that 
containing the incident and scattered beams. Then, 

Ei = ai exp(-i$i) exp(i(mt-k^. (2.21) 

E r = ai exp(-i5 r ) exp(i(ut"kz), 

6 = — 6j- 

As an alternative to requiring that the amplitudes A^ and A r , and the 
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phase difference $ be known, the state of polarization may be represented 

by the four Stokes parameters Introduced by Sir George Stokes In 1852! 

These have the advantage of all having the same dimension, that of an 

irradiance. The four parameters are 

II - E1E1* - 01* ( 2 . 22 ) 

I r - ErE r * “ a r * 

U » 2 Re(EiE r *) ** 2 aia r cos6 
V “ 2 Im(EiE r *) “ 2 aia r sln<S 

where the asterisk denotes the complex conjugate has been taken. 

Referring to Figure 2.3, the state of polarization can be 
represented by an ellipse, which In turn Is described by the Stokes 
parameters. Let x be the angle between the direction of the major axis 
and the 1 direction. Knowing x is equivalent to knowing the plane of 
polarization, or that plane through Che direction of propagation and ray 
containing the maximum electric field vector. Also, let the ellipticity 
be represented by the angle @ whose tangent is the ratio of the lengths 
of the major and minor axes. It can be shown, as in Chandrasekhar (1950), 
that 

I - Ii + I r (2.23) 

Q = II - Ir 

tan2x 3 U/Q 
sin28 - V/I 

Therefore, the parameters !]_, I r , U, and V represent the irradiances in 
two perpendicular directions within a plane transverse to the direction of 
propagation, the plane of polarization, and the ellipticity of the 
electromagnetic wave. With these, all quantities relevant to the 
description of the state of polarization are determined. In addition, the 
percentage polarization is given as 

P = /Q*+U 2 -rV 2 /l (2.24) 
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Figure 2.3 Representation- of elliptical polarization 
(from Llou, 1980). 




Figure 2.4 Dipole scattering. 
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For unpolarized light U“V=0, as the time average over sin6 and cos6 are 
zero, and Ii*I r . The polarization P equals zero, as expected. One 
example of unpolarized light is that which is initially from the sun, 
although the light becomes partially polarized after scattering within the 
atmosphere. Conversely, for a completely polarized beam I 1 =Q 1 +U I +V I , and 
P-1. 


The expressions 


% - Si Eol (2.25) 

E r = Sj Eox 

are next utilized to determine the four Stokes parameters of the 
scattered field. Si and S r relate the magnitude of the scattered fields 
Eiand E r to that of the incident fields Eol anc * Eor> They are functions 
of the angle between the incident and scattered directions of propagation, 
and, as will be shown later, differ amongst the Rayleigh and Mie 
particles. By substituting Equation (2.25) into (2.22), the four Stokes 
parameters are determined 

Tl n Si Si* I.i (2.26) 

Ij- = S r S r lot* 

u = u 0 Re(SiSr*) + V 0 Im(SiS r *) 

V = Uo Im(SiS r *) + Vo Re(SiS r *) 

In much the same way as Equations (2.16a,b) represents the 
transfer of radiant flux within the atmosphere, they likewise can 
represent the transformation of the Stokes parameters as the beam they 
represent undergoes scattering within the atmosphere. A few 
modifications are, however, required. First, from Equation (2.26) it can 
be shown that not all the scattered Stokes parameters are independent of 
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each other. To account for this, the radiative transfer equation is 
rewritten in matrix form. That is, 


L pC T f * ± w, $) 


L p (T lt ±u,+) e-At/y + 


T n+Ar 

JpqCf'jiy, $) e“ iT, /p dt'/y 
T n 


(2.27a) 


and 


J pqC T f * ± M, 4>) 


p pq(y , »4' , ;±y»<l') Lq(T , »M , »*') doj' 
4 tt 


(2.27b) 


+ PpqC-U 0 , 4> 0 ;±P» 4>) E 8 q e-T'/Po . 

Here p and q are related to one of the four Stokes parameters, and Ppq is 
a 4x4 matrix. Thus, the p c ^ component of radiance is determined by 
summing the source function over the four incoming Stokes components, i.e. 


q 13 ! ,2,3,4. 


The matrix Ppq cannot be written directly from Equation (2.26). 
As the equation of radiative transfer traces components relative to a 
vertical plane within the atmosphere, not the scattering plane, a 
coordinate transformation must be performed. The most general such 
scattering phase matrix has been given by Sekera (1955). It takes the 
form 


Ppq = 


A i i A n* 

A 2 iAj i 

2 Re(A ll A z j*) 
2 Im(A lt A zl *) 


A A * 

A, 

A A * 

A 22 A 22 . 

2 Re(A l2 A 22 *) 
2 Im(A 12 A 22 *) 


Re(Aj jA 12 *) -Im(A xl A l2 *) 

Re(A 2l A 22 *) -Im(A 2l A 22 *) (2.28) 

Re(A, 2 *A Z1 +A, ,A 2I ) — Im(A 21 Aj 2 ^" A i * A 22 

Ira(Aj x A 22 -A 12 A 21 ) Re(Aj jA 22 — A 2 |A 12 ) 


The quantities Apq are given by 

An =11 cosA(|) + T 2 cosip (2.29) 

A l2 “ (y'Ti + yT z ) sinA<j) 

A 2 1 = (yTj + y'T z ) sinA$ 

A 22 “ T l COS<|> + T 2 COS A c|l . 
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Here p and p' are directional cosines of the incident and scattered beams 
(measured, as before, from the local vertical). The angle ia defined 

as the difference between the azimuth angles of the incident and 

scattered beams. Furthermore, 

cos 41 ** (l~vi 1 ) l / J (l-n ,z ) 1 /* + y|i' cos &,{> (2.30) 

T t - (S X - XS r )/(l-X 4 ) 

T* » (S r - XSi)/(l-X*) 
x « cose » nu' + (i-p*) l /*(l-V*) 1 / 1 cos a<j> 

and Si, S r are the proportionality constants defined in Equation (2.25). 

It is to be noted that all the functions within Equations (2.28) through 

(2.30) are defined with respect to the angles 0 , 0 ’, and 

Rayleigh Scattering by Molecules 
Both molecules, whose size are on the order of 10“^ pm, and 
aerosols, ranging from 0.01 to 10 Pm, are responsible for scattering 
within the atmosphere. Molecular scattering in the visible and near ir, 
where 2ur<<X, can be characterized by a simple scattering law due to Lord 
Rayleigh (J.W. Strutt, third Baron of Rayleigh). In 1872 be derived the 
scattering law, which now bears bis name, using the elastic-solid ether 
theory. He predicted that scattering varies inversely as the fourth power 
of the wavelength, and so explained the blue color of the sky. In 1899 
Rayleigh revised his derivation to use the electromagnetic theories of 
Maxwell and Hertz, Thus, the dependence of scattering on refractive 
index was determined. The scattering law has since undergone one slight 
revision, to account for molecular anisotropy. This was done in the 
1920's, shortly after some scattering experiments made by Rayleigh's son 
demonstrated the need for this modification. A complete development of 
the Rayleigh scattering law is given in texts such as McCartney (1976). 
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Highlights of Its development are given here 


Dipole Scattering 

To begin with, the mechanical oscillator model of the atom is 
used. A binding force is characterized by a spring which induces a linear 
restoring force to the electron as it is displaced. Such a displacement 
occurs when a molecule is subject to an applied electric field, E 0 , An 
induced dipole moment p«ex ia created, where e is the charge on an 
electron, and x is the displacement. This electric dipole oscillates 
synchronously with the field, and in turn produces the scattered wave. 
The new field is proportional to (l)tbe acceleration of the electron, 
(2)sln6, where e is the angle between the dipole moment and direction of 
observation, and it is inversely proportional to R, the distance from the 
dipole. It has an amplitude 

cu 2 p 0 sine sinw(t-R/c) 


4ire 0 C “ R 


(2.31) 


Because of the sin0 dependance, the dipole cannot radiate along the axis 
of the dipole. The maximum dipole moment p„ is found by solving the 
equation of motion for the maximum electron displacement: 


e 2 E. 


e x 


(n 2 -l) 3 e o E 0 


0 m(<o 0 2 -*ii) 2 ) 


(2.32) 


(n 2 +l) N * 

Here, o}„ is the resonant frequency of oscillation, equal to (k/m) 1 /* where 
k is the restoring force on the electron. 

The latter equality within Equation (2.32) utilizes the Lorenz- 
Lorentz expression to substitute for the molecular parameters. Now n, 
the refractive index of the gas in bulk form, and N, the number of dipole 
oscillators per unit volume, are used. The refractive index of air 
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molecules considered here are found Co be nearly 1, real, and vary as a 
funcCion of wavelength. This wavelength dependence is given by Edlen 
(1953) as 


(n-1) x 10* - 6432.8 + • 2 - ' 9 - 4 ^ 810 + 25*540 

1 Afi-l -* A1-1-* 


(2.33) 


146-X -2 41-r 

For example, 0 “ 1 . 000293 at X-0.55 pm. Both the Lorenz-Lorentz and Edlen 
expressions are derived in many discussions on the dispersion of 
electromagnetic waves, as in Liou (1980). 

The Rayleigh expressions assume that scatterers have resonant 
frequencies far above the visible and infrared spectral regions. Thus 
they are pure scatterers, and absorb no energy. Such an assumption is 
valid for nitrogen and oxygen molecules, which are responsible for 99 % of 
molecular scattering. There are, however, molecules that do have an 
imaginary component to their refractive index at those wavelengths of 
interest (i.e., they have resonant frequencies near those frequencies 
corresponding to visible light). The effects of scattering from these 
species can be overlooked without loss of accuracy, as they compose such 
a small fraction of the atmospheric gases. Ozone and water vapor are two 
such absorbers. (The columnar amount of ozone is typically on the order 
of 0.35 cm-atm. This implies that there will be only 0.35 cm of ozone 
within a 1 cm 2 atmospheric column of air, in which there are several 
kilometers of atmospheric scatterers.) 

The irradiance produced at a distant point R from the dipole is 
given by the Poyning vector "S’, 


S - c e 0 <E 2 > . (2.34) 

The mean of E 2 is found by substituting a factor of 1/2 for sin 2 w(t-R/c), 
and using Equations (2.31) and (2.32) for the electric field strength. To 
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remove the dependence of scattering on distance R, the intensity, I 
(Watts/sr), is computed instead. The intensity at distance R is found by 
multiplying the irrsdiance S by R 2 (since I=d<i>/dw“d<|>/dA x dA/dw, and dA"R z 
dw). Hence, 

_ w 2 e 0 c sin 2 © (n 2 -l ) 2 E 0 2 /0 oc , 

( 0 ) 2 N 2 \ H * (2.35) 

In addition to the previous equations, the substitution w=2irc/X and 

(n+2)*“9 (since n^l) have been mode in writing Equation (2.35). 

Cross Section 

The scattering cross section of a gas molecule is defined as that 
cross section of an incident wave, acted on by the molecule, having an 
area such that the irradlance flowing across it is equal to the total 
irradiance scattered in all directions. Thus, 


/ 

I( 0 ') du' 
4 7T 

c"c, E //2 


(2.36) 


Using Equations (2.35), du'=2r sin©' d0', and 


sin*©' d ©'=4/3, 

, 0 


the cross 


section is obtained. To this, the correction factor (6+36)/(6-76) must be 
added. This is done to account for molecular anisotropy, which prevents 
the dipole moment from aligning itself exactly with the electric vector 
of the primary wave. Thus, 


8 ir 3 (n 2 ~l ) 2 6+3 6 
ff Ray = 3 N i V 6-7<S * 

Gucker and Basu (1953) have determined that 6=0.035. 


(2.37) 
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Rayleigh Optical Depth 

The volume scattering coefficient for molecules, $Ray* gives the 
fractional amount of flux scattered in all directions, for a unit volume 
of gas. As the scattered field from a collection of dipoles add 
incoherently, the angular coefficient for a unit volume is Just N 
(molecules/volume) times the cross section given by (2.37), or $Ray*crRay 
N. (Also, the mass extinction coefficient is found to be kRay*°Ray 
N/p"<jRay / m , or cross section per unit mass, where m is the moss of the 
molecule). Using the definition of optical depth, Equation (2.9), the 
Rayleigh component of optical depth is determined 


tgay » a Ray 


f N(z) dz. 


(2.38) 


Jz 

Model values of the molecular number density as a function of altitude 
can be found in the U.S. Standard Atmosphere-1962 (see, for example, 
Valley (1965), or Elterman (1968)), and are given here in Table 3.2. 

The tabulated values of mass density, p, or number density, N, 
refer to air at sea-level temperature and pressure. It Is desirable to 
compute the scattering coefficients at nonstandard values of temperature, 
pressure, and altitude. This is done using the equation of state for an 
ideal gas (P»pRT, P being atmospheric pressure, R the universal gas 
constant, and T the temperature on the Kelvin scale). Thus, 


P T 0 


p o P # T * 


(2.39) 


where p 0 , P 0 , and T 0 are defined at standard atmosphere conditions. 

In using the Herman code to model the atmosphere, tRay is 
determined using measured values of atmospheric pressure. At ground 
level, . 
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V 


T Ray 


8ir J (nM)* 

3 ~X k N 1 


Nr 


10 1 ' 



(2.40) 


where 

n - * >frai:tive index ns given by Equation (2.33) 

X « „i ve length in Mm 

Ng * molecular number density at sea level for a standard 
atmosphere 

- 2.547 x 10*' cm**’ 

N c - columnar number density 

- 2.154 x 10* ‘ cm" 1 
d - 0.035 

P 0 “ 1013.25 mbar, or 29.9? in Hg 

P » measured atmospheric pressure, same units as Po. 

Using this formulism at A-0,55 Mm, for example, TR a y*.098, 


Phase Function 

The angular dependence on scattering is expressed in terms of the 
phase function P(0). This function is defined as the ratio of the radiance 
into a given direction, to the average radiance in all directions. Thus, 
the integral of the phase function must be normalised to unity, as there 
is no absorption by Rayleigh molecules, and 


/ 

P(0) dw* » <o 0 - 1 . (2.41) 

4ir 

To derive the phase function for the scattering of unpolarized 
light by Rayleigh particles, the incident electric field vector is 
decomposed into two orthogonal components. As before, let E^ and E r 
represent those scattered components parallel^ and perpendicular^ to a 
reference plane, and let E 0 1 and Eor be the corresponding incident 
components. The reference plane is taken as that containing the incident, 
and scattered waves, and the scattered wave is deviated from the incident 
wave by an angle 8. For each of these two components (i=l or r), the 
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scattered radiance is found using 


L(0i) ■ 1(0) (6+34)/(6-76) N ds - 3/8* oj^y N ds sin 1 ©! S*i (2.42) 
In arriving at this expression, it Is noted that as 1(0) gives the 
intensity scattered from a single molecule, I(0)N do gives the radiance 
scattered from a volume of gas. After accounting for anisotropy, the 
intensity is expressed in terms of OR a y by using (2.37). The incident 
irradiance cc 4 E 0 */2 is then expressed as S s !» 

The angles 0i and 0 r can readily be expressed in terras of Che 
scattering angle, 9. With reference to Figure 2.4, it is shown that 
01 - ir/2~Q, and 0 r “ir/2. Hence, the total scattered radiance is given as 

L " L r + Li “ 3/8 it o N ds L #r + 3/8w a N ds 008*0 S, r . (2.43) 

But, as the incoming field is unpolarized, L Br “Loi“L 9 /2. Equation (2,43) 
becomes 

L * o N ds 3/ 1 6 it (1 + cos 2 0) s 0 . (2,44) 

Removing Che angular dependence and multiplying by a scaling factor to 
satisfy (2.41), the phase function for Rayleigh scattering of unpolarized 
light is found to be 

P(9) - 3/ 1 6 tt (1 + 008*9 ) . (2.45) 

This expression is the Rayleigh component of the phase function used 
within Equation (2.16). It is thus an important parameter in the 
calculations of the transfer of radiant flux within the & Vmosphere. 

Because the perpendicular and parallel components are not 
scattered equally, the resulting radiance will be partially polarized. 
Although the scattered perpendicular component is independent of Che 
angle 0, the parallel component follows a cos 2 9 dependance. Thus, if the 
observation direction is at 90° to the incoming beam, the scattered light 
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will be completely polarized# The scattered energy is symmetric about 
the incident beam, and equal amount of energy are sent into the forward 
and backward hemispheres. If there were only single scattering within 
the atmosphere, and the atmosphere was composed purely of Rayleigh 
particles, the skylight everywhere from a 90° angle from the earth-sun 
line would be completely polarized. This perfectly polarized light is 
never observed in practice, as the scattering from aerosols, the reflected 
light from the surface, and the anisotropby of air molecules themselves 
cannot be neglected. 

If these polarization changes are to be traced through the 

atmosphere, the matrix form of the phase function is required. For 

Rayleigh scattering, this becomes (2,46) 

cos 2 i|j p 2 sin 1 i(j) pcoe vp sinA$ 0 

Ppq “ 3/8 tt u' a sirt 2 A 4* cos 2 4«t> -p'sinA^ cosA4> 0 

-2u'co8'{) si , nA‘{> 2psinA4> cosA<}> -jiji'sin 2 A4>+cosil> co8A$ 0 

0 0 0 costy cosA«|»+yiJ , 8in 2 A<j) 

Mie Scattering 

To describe scattering by particles of arbitrary size the 
equations developed by Mie (1908) are universally used. In developing 
this theory it was necessary to make the simplifying assumption that the 
scattering particles were isotropic spheres. Even so, the derivation is 
complex, using Maxwell's equations, a boundary value analysis, and 
expansion of the emerging wave in terms of a series of Bessel and 
Legendre polynomials. The equations can be approximated by the first 
term of the Mie series for small particles. For this case, however, 

Rayleigh theory yields an equivalent result with significantly fewer 
computations. Thus, the term Mie scattering is loosely used to refer to 
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the scattering by larger particles which do not lie within the Rayleigh 
regime. 

A complete development of Mle theory Is given by Stratton (1941) 
and van de Hulst (1957). The scattered light Is again found by breaking 
the Incident beam Into components perpendlcula, and parallel to the 
scattering plane. The scattered Intensities I r (0) and li(0) are 
proportional to the functions 

09 

i r * |2ir/X S r | 2 » | ^ C a n 1T n +lj n T n) | 2 (2.47a) 

n«l 

90 

*1 - I 2tt/ X Si | ; - | JT < a n T t>+1>ii*n> I* (2.47b) 

n»l * 

Each function Is found as the sum of an Infinite series. Defining the size 
parameter as a=2irr/X where r Is the the radius of the particle, It Is 
found that the number of terms required for convergence Is somewhat 
greater than a, for a>l. The amplitudes of the nth electric partial wave 
and the nth magnetic wave are given by the complex coefficients an and 
bp. These are 

jn(a)J' ~ jn(a)[ma jn( m a)]' 

a n = ~ — — (2.48a) 

jn(ma)[a bn^JCa)!' - b n ( 2 )(a)[ma jn(ma)]' 

jnCa)[ma inC^a)!' ” m* jnC m a)Eci j nC a ) ] * 
b n = . (2.48b) 

bfl( 2 )(a)[nio jn(ma)]' - m 2 jn(ma)[a b n ^Ha)]' 

With air as the Incident medium, the parameter m= , n re (l-ni !n ^ related to 
both the real and Imaginary components of the refractive Index within the 
sphere. Spherical Bessel and Hankel functions are denoted by jn and bn 
respectively, and primes denote derivatives with respect to the indicated 
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arguments. Thus the coefficients a n and b n are determined from the 
particle characteristics, but are independent of the scattering angle 9. 
This latter dependance is expressed through the functions ir n and x n and 
involve the first and second derivatives of Legendre polynomials: 

d(P n (cose)) 


irn(cose) 


d(cosg) 


(2.49a) 


Tn(cose) a cos 9 ir n (cos9) - sin 4 9 (2.49b) 

When the particle is illuminated by plane-polarized light, the 
intensity of the scattered light is given by 

1(9) a E 0 (i r sin 2 \]j + ii cos 2 ^) . (2.50) 

Here E 0 is the irradiance of the incoming beam, ^ is the angle of the 
electric vector from the scattering plane, and i r and li are as defined in 
(2.47). For a particle illuminated with a wave whose electric vector is 
perpendicular to the plane of observation, f=90° and the scattered beam 
is polarized in the perpendicular direction. Conversely, an incident beam 
described by ^=>0 is polarized parallel to the scattered plane, as is the 
scattered beam. For illumination by an unpolarized beam, the scattered 
intensity is given by 


1(e) - E or Ir + Ed il (2.51) 

' E * & ( f r + il) 

where E 5T =E 0 i=E fl /2 . 

The angular distribution of the scattered field is depicted in 
Figure 2.5. Here the solid lines refer to scattering from a perpendicular 
component of tbe electric vector, and the dashed lines represent 
scattering from a parallel component. For a<0.1 tbe distribution is 
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Figure 2.5 Mie scatteri 
(from Grama, 1978). 
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Identical to that predicted from Rayleigh theory. There is a cos 2 9 
dependance in the scattered parallel component but no angular variation 
in the perpendicular component. As a increases (or particle size for a 
given wavelength), a larger portion of the energy is scattered into the 
forward direction. If the particle size approaches the wavelength of 
light, side lobes begin to appear. The frequency of this structure 
increases with a and the width decreases. 

The cross section of a Mie scatterer can now be defined. Unlike 
scattering from a Rayleigh particle, some energy is lo^t due to absorption 
as a beam impinges upon a Mie scatterer. The cross section muBt 

include the effects of this loss.' Defining a 8C as the component which 
accounts for the energy scattered into all directions, and o a t, 8 the 
component which accounts for absorption, we have 

a Mie = °sc + ff abs • (2.52) 

Using Equation (2.50) and assuming unpolarlzed illumination, o sc is 
computed from 


< T sc “ 


f 


1(9) dm/E 0 = A 2 / 87r 2 | (i r + ij)sino de d$ 


(2.53) 


" * 2 / 2 * JT (2o+l)(ia n 1 2 + lb n 5 2 ) . 
n=l 

The total cross section can likewise be expressed in terms of the Mie 
coefficients, 


oMie - X 2 / 2 tt JT(2n+l) Re(a n + b n ) . (2.54) 

n=l 

In the above the expansions in terms of an and bn do not easily follow. 
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Reference is given to van de Hulst (1957, section 9.32) for more details. 

It is noted that the above cross sections are defined for a 
particle of fixed radius r. Absorption and scattering within the 
atmosphere are, however, processes which depend on the cumulative effects 
of many particles within a large size range. This distribution is 
expressed in terms of a size distribution n(r). It is the normalized 
number of particles per unit interval of radius per unit volume, hence 

/• 

n(r) dr * 1 . (2.55) 

0 

To determine the properties of light scattered from a polydispersion 
(collection of particles of different radii), the functions i r and ij 
within (2.53) are integrated over the size distribution. The scattered 
energy from such a distribution of particles is very different from that 
depicted above. The most obvious difference is that the scattered 
distribution is a much smoother function of wavelength. A few examples 
of this are given in Figure 2.6. To compute these curves a log-normal 
particle distribution was assumed. A mean radius of r m =l pm, standard 
deviation 2 pm, wavelength X ct 0.633 pm, and real refractive index n re = , 1.525 
were assumed. Curves (a) and (b) give the the results for a parallel and 
perpendicular incident electric vector, respectively. In curves (c) and 
(d) the molecular scattering contributions have been added. Each example 
has been computed at several values of the imaginary component of 
refractive index. As n-j m increases the light scattered into angles 
greater than 9=15° decreases. The most significant result of increasing 
the imaginary refractive index, however, is the increase in absorption. 

This change can be expressed through the parameters a a bs or id 0 , the 
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Figure 2.6 Mie scat 
(from Grams, 1978). 





t 


single scatter albedo. 

By Integrating the cross section over the size distribution 
function, the optical depth of the atmosphere can be determined. Defining 
N(z) as the totax number of particles per unit volume at altitude z, 


T Mie 


N(z) o m ie( r ) °( r ) dr dz * 


(2.56) 


Note that the size distribution is taken as constant with respect to 
altitude. This is usually assumed the case, for lack of better data. 
More will be said about the radial size distribution function n(r) and the 
vertical distribution N(z) in the sections to follow. 

By integrating the cross section over the size distribution, the 
phase function for Mie scattering can also be found. To most readily see 
this, let us define the angular scattering cross section agc(9) as the 
cross section of the incident wave acted on by the particle, having an 
area such that the irradiance flowing across it is equal to the intensity 
scattered into angle 0. The cross section o sc defined earlier is equal to 
the angular cross section integrated over all outgoing angles. With this, 
the phase function is defined as 


' 

a sc( 0 >O n(r) dr 

P(0) - . (2.57) 

| <mie(r) n(r) dr 

From this definition it is apparent that the integral of the phase 
function over all solid angles will not necessarily be equal to one. It 
will be equal to ti> 0 , the single scatter albedo, and equal to one only if 
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there Is no absorption of energy by the particle (In such a case 
0Mle"°sc)* The greater the Imaginary component of refractive Index, the 
smaller u> fl will be, bence a smaller fraction of energy will be scattered, 
a greater fraction absorbed. 

For Illumination by unpolarized light 


P(e) 


*V8tt j 


n(r)(l r + ii) dr 


[ aMle(r) n(r) dr 


(2.58) 


To run that version of the Herman Code which accounts for polarization, 
the phase function must be written In matrix form. This Is done by using 
S r and Sl> as defined in (2.47), within Equations (2.28)~(2.30). 

Another parameter, closely related to the cross section, that is 
commonly referred to in the literature is the efficiency factor Q, defined 
as the cross section of a particle divided by the geometric cross- 
sectional area of that particle, irr 2 . If the scattering efficiency factor 
is ploted versus the size paramter a, Q gc obtains a maximum value of 2 
and converges in an oscillatory fashion to a value of one for high a. 
This implies that the particle can, at times, interact with an incident 
wavefront greater than its own geometric area. This is explained through 
diffraction effects, in which diffracted flux is directed into a small 
angle centered about the forward direction of the incident flux. 


Dave Code 

To compute the Mie parameters discussed above, a Fortran computer 
program written by Dave (1969) is used. This program is incorporated 
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into one of the subroutines within the Herman Code. Two slmiliar Dove 
codes exist, one using an upward recurrence relationship (starting with a 
value of a 9 (ma), successively higher values are computed), and one using a 
downward recurrence relationship. In Che code which uses an upward 
recurrence algorithm, any error In the first term will propagate and for 
large enough ot the results oscillate wildly around the correct value. 
For this reason the downward recurrence routine is preferred. It does, 
however, more storage, and required 10-20% more run time. Both codes 
require double precision arithmetic, and output results accurate to 6 
significant figures (with the one exception mentioned above, where 
oscillations occur). 
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